% Compute the global solution using a paramerized expectations algorithm
% Dynamic bank run model
% Bank leverage choice only
% Capital depreciation shock

close all
warning off

load('rbc')

%% Generate random shocks

options = optimset;
T = 10000;               % simulation periods
t_drop = 1000;

rng(1)                   % same random shock 
shocks = normrnd(0,1,T+1,1);  % random shock
%Exi = zeros(T+1);  % stochastic steady state

kv   = zeros(T+1,1);
rbv  = zeros(T+1,1);
nv   = zeros(T+1,1);
xiv  = zeros(T+1,1);
hv   = zeros(T+1,1);
cv   = zeros(T+1,1);
yv   = zeros(T+1,1);
iv   = zeros(T+1,1);
lv   = zeros(T+1,1);
pibv = zeros(T+1,1);
runv = zeros(T+1,1);
rhseev = zeros(T+1,1);
rkhat  = zeros(T+1,1);
rkhats = zeros(T+1,1);
prv  = zeros(T+1,1);
av   = zeros(T+1,1);
erk  = zeros(T+1,1);
x    = zeros(T+1,npol);
x_nonlog    = zeros(T+1,npol);
x_allnonlog = zeros(T+1,npol);

%% Simulation

kv(1)  = k_ss;
av(1)  = a_ss;

    for t=2:T+1
        av(t) = rho_a*av(t-1) + sigma_a*shocks(t);
        hv(t)  = ((1-alpha)/psi*exp(av(t))*kv(t-1)^(alpha))^(1/(alpha+1/nu));
        yv(t)  = exp(av(t)) * kv(t-1)^(alpha)*hv(t)^(1-alpha);

        x0 = [log(kv(t-1)) av(t)].^expmx;
        x(t,:)  = prod(x0,2)';
        
        rhseev(t) = exp(x(t,:)*solved_eta1_rhsee);
        cv(t) = rhseev(t)^(-1) + psi*hv(t)^(1+1/nu)/(1+1/nu);
        iv(t) = yv(t) - cv(t);
        kv(t) = (1-delta)*kv(t-1) + iv(t);

        rbv(t) = 1 - delta + alpha*yv(t)/kv(t-1);

    end
    % End of T period simulation

c_alltime_avg  = mean(cv(t_drop+1:T));
y_alltime_avg  = mean(yv(t_drop+1:T));
h_alltime_avg  = mean(hv(t_drop+1:T));
k_alltime_avg  = mean(kv(t_drop+1:T));
i_alltime_avg  = mean(iv(t_drop+1:T));

c_alltime_std  = std(cv(t_drop+1:T));
y_alltime_std  = std(yv(t_drop+1:T));
h_alltime_std  = std(hv(t_drop+1:T));
k_alltime_std  = std(kv(t_drop+1:T));
i_alltime_std  = std(iv(t_drop+1:T));

return
